##
## Distance to roads
##

require(Rcpp)
output_version <- '2020-09-30'
options(rasterTmpDir='C:/temp/')
ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones)


##
## Roads
##
road_lk_chad_shp = paste0(wkdir,'/local/gis_data/018_transportation/roads_panel')
roads <- st_read(dirname(road_lk_chad_shp),layer=basename(road_lk_chad_shp))
road_cat <- sort(unique(roads$R2008))

ntl.zones.pts$dist2road0<-NA
ntl.zones.pts$dist2road1<-NA
ntl.zones.pts$dist2road3<-NA
ntl.zones.pts$dist2road5<-NA
ntl.zones.pts$dist2road8<-NA

# CMR
ntl.zones.pts$dist2road0[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], subset(roads,R2008==0 & COUNTRY=="Cameroon"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road1[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], subset(roads,R2008==1 & COUNTRY=="Cameroon"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road3[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], subset(roads,R2008==3 & COUNTRY=="Cameroon"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road5[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], subset(roads,R2008==5 & COUNTRY=="Cameroon"), k = 1, returnDist = T)$dist) / 1000

# NER
ntl.zones.pts$dist2road1[ntl.zones.pts$GID_0=="NER"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NER",], subset(roads,R2008==1 & COUNTRY=="Niger"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road3[ntl.zones.pts$GID_0=="NER"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NER",], subset(roads,R2008==3 & COUNTRY=="Niger"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road5[ntl.zones.pts$GID_0=="NER"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NER",], subset(roads,R2008==5 & COUNTRY=="Niger"), k = 1, returnDist = T)$dist) / 1000

# NGA
ntl.zones.pts$dist2road1[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], subset(roads,R2008==1 & COUNTRY=="Nigeria"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road3[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], subset(roads,R2008==3 & COUNTRY=="Nigeria"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road5[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], subset(roads,R2008==5 & COUNTRY=="Nigeria"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road8[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], subset(roads,R2008==8 & COUNTRY=="Nigeria"), k = 1, returnDist = T)$dist) / 1000

# TCD
ntl.zones.pts$dist2road1[ntl.zones.pts$GID_0=="TCD"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="TCD",], subset(roads,R2008==1 & COUNTRY=="Chad"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road3[ntl.zones.pts$GID_0=="TCD"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="TCD",], subset(roads,R2008==3 & COUNTRY=="Chad"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts$dist2road5[ntl.zones.pts$GID_0=="TCD"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="TCD",], subset(roads,R2008==5 & COUNTRY=="Chad"), k = 1, returnDist = T)$dist) / 1000

save.dta13(st_drop_geometry(ntl.zones.pts), paste0(wkdir,"/proc_data/DIST2ROAD08_",output_version,".dta"))